Phase transition approach to detecting singularities 

of PDEs 



Panagiotis Stinis 
Department of Mathematics 
University of Minnesota 
Minneapolis, MN 55455 

Abstract 

We present a mesh refinement algorithm for detecting singularities 
of time-dependent partial differential equations. The algorithm is in- 
spired by renormalization constructions used in statistical mechanics 
to evaluate the properties of a system near a critical point, i.e. a phase 
transition. The main idea behind the algorithm is to treat the occur- 
rence of singularities of time-dependent partial differential equations 
as phase transitions. 

The algorithm assumes the knowledge of an accurate reduced model. 
In particular, we need only assume that we know the functional form 
of the reduced model, i.e. the terms appearing in the reduced model, 
but not necessarily their coefficients. We provide a way of computing 
the necessary coefficients on the fly as needed. 

We show how the mesh refinement algorithm can be used to calcu- 
late the blow-up rate as we approach the singularity. This calculation 
can be done in three different ways: i) the direct approach where one 
monitors the blowing-up quantity as it approaches the singularity and 
uses the data to calculate the blow-up rate ; ii) the "phase transition" 
approach (a la Wilson) where one treats the singularity as a fixed 
point of the renormalization flow equation and proceeds to compute 
the blow-up rate via an analysis in the vicinity of the fixed point and 
iii) the "scaling" approach (a la Widom-Kadanoff) where one postu- 
lates the existence of scaling laws for different quantities close to the 
singularity, computes the associated exponents and then uses them to 
estimate the blow-up rate. Our algorithm allows a unified presentation 
of these three approaches. 

The inviscid Burgers and the supercritical focusing Schrodinger 
equations are used as instructive examples to illustrate the construc- 
tions. 



Introduction 



The problem of how to construct mesh refinement methods and how to ap- 
proach more efficiently possible singularities of partial differential equations 
has attracted considerable attention (see e.g. [H HI [31 [H HE]). At the 
same time, the problem of constructing dimensionally reduced models for 
large systems of ordinary differential equations (this covers the case of partial 
differential equations after discretization or series expansion of the solution) 
has also received considerable attention e.g. see the review papers |12l lllj . 
The construction of an accurate reduced model has advantages beyond the 
obvious one of predicting the correct behavior for a reduced set of variables. 

We present here an algorithm that is based on dimensional reduction 
and which can be used to perform mesh refinement and investigate possi- 
bly singular solutions of partial differential equations (see also [18]). The 
algorithm is inspired by constructions used in statistical mechanics to eval- 
uate the properties of a system near a critical point [131 [5] (a critical point 
is a value for the controlling parameter of a system at which the behavior 
of the system changes abruptly). The idea underlying the computation of 
the properties at criticality, is that while the form of the reduced system 
equations is important, one can extract even more information by looking 
at how the form of the reduced system is related to the form of the original 
(full dimensional) system [2Ql [21] . We extend this idea to the study of (pos- 
sibly) singular solutions of partial differential equations by treating time as 
the controlling parameter and the instant of occurrence of a singularity as 
a critical value for the parameter, i.e. a critical point. 

Our approach has two objectives: i) it provides a way of accurately mon- 
itoring the progress of a simulation towards underresolution, thus providing 
us as a byproduct with the time of occurrence of the possible singularity; 
ii) it allows the formulation of a mesh refinement scheme which is able to 
reach the time of interesting dynamics of the equation much more efficiently 
compared to an algorithm that simply starts with the maximum available 
resolution. 

The mesh refinement algorithm can be used to calculate the blow-up 
rate as we approach the singularity. This calculation can be done in three 
different ways: i) the direct approach where one monitors the blowing-up 
quantity as it approaches the singularity and uses the data to calculate 
the blow-up rate ; ii) the "phase transition" approach (a la Wilson) [13] 
where one treats the singularity as a fixed point of the renormalization flow 
equation and proceeds to compute the blow-up rate via an analysis in the 
vicinity of the fixed point and iii) the "scaling" approach (a la Widom- 



2 



Kadanoff ) [5\ where one postulates the existence of scaling laws for different 
quantities close to the singularity, computes the associated exponents and 
then uses them to estimate the blow-up rate. Our algorithm allows a unified 
presentation of these three approaches. 

The task of investigating numerically the appearance of a singularity is 
subtle. Clearly, since all calculations are performed with finite resolution 
and a singularity involves an infinity of active scales we can only come as 
close to the singularity as our resolution will allow. After that point, either 
we stop our calculations and conclude that a singularity may be present 
close to the time instant we stopped or we switch, if available, to a model 
that drains energy at the correct rate out of the set of resolved variables. We 
emphasize here that, up to some time, the evolution towards a near-singular 
solution can be identical to the evolution towards a singular solution. If we 
do not have enough resolution to go beyond the time instant after which the 
two evolutions start deviating, we cannot claim with certainty the presence 
of a singularity. In other words, given adequate resolution we can eliminate 
the possibility of a singularity. But, it may be very hard, to prove through 
a finite calculation, that a singularity exists (we come back to these points 
in Sections [2] and \3$) . 

The paper is organized as follows. In Section [1] we present the ideas 
behind the construction of the algorithm. In Section [2] we present the mesh 
refinement algorithm. In Section [3] we provide numerical results for the 
inviscid Burgers equation. In Section [J] we provide numerical results for the 
supercritical focusing Schrddinger equation. Section [5] shows how one can 
use the mesh refinement algorithm to compute the blow-up rate as a critical 
exponent, i.e. using solely properties of a renormalization (coarse-graining) 
process in the vicinity of the singularity. Section [6] contains a discussion of 
the results and some directions for future work. 

1 The main construction 

Suppose that we are interested in the possible development of singularities 
in the solution v(x,t) of a partial differential equation (PDE) 

v t + H(t,x,v,v x , ...) = 

where H is a, in general nonlinear, operator and x £ D C M. d (the construc- 
tions extend readily to the case of systems of partial differential equations). 
After spatial discretization or expansion of the solution in series, the PDE 
transforms into a system of ordinary differential equations (ODEs). For 
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simplicity we restrict ourselves to the case of periodic boundary conditions, 
so that a Fourier expansion of the solution leads to system of ODEs for the 
Fourier coefficients. To simulate the system for the Fourier coefficients we 
need to truncate at some point the Fourier expansion. Let FUG denote the 
set of Fourier modes retained in the series, where we have split the Fourier 
modes in two sets, F and G. We call the modes in F resolved and the modes 
in G unresolved. One can construct, in principle, an exact reduced model 
for the modes in F e.g. through the Mori-Zwanzig formalism [9] (we do not 
deal here with the complications of constructing a good reduced model) . 

The main idea behind the algorithm is that the evolution of moments of 
the reduced set of modes, for example L norms of the modes in F, should 
be the same whether computed from the full or the reduced system. This 
is a generalization to time-dependent systems of the principle used in the 
theory of equilibrium phase transitions to compute the critical exponents 
[131 [T7] . The idea underlying the computation of the critical exponents is 
that while the form of the reduced system equations is important, one can 
extract even more information by looking at how the form of the reduced 
system is related to the form of the original (full dimensional) system. We 
extend this idea to the study of (possibly) singular solutions of partial dif- 
ferential equations by treating time as the controlling parameter and the 
instant of occurrence of a singularity as a critical value for the parameter, 
i.e. a critical point. We caution the reader that even though our motiva- 
tion for the present construction came from the theory of equilibrium phase 
transitions, we do not advocate that a singularity is a phase transition in 
the conventional sense. It can be thought of as a transition from a strong 
solution to an appropriately defined weak solution but one does not have 
to push the analogy further. We want to point here that the problem we 
are addressing is different from the subject known as dynamic critical phe- 
nomena (see Ch. 8 in |13]). There, one is interested in the computation 
of time-dependent quantities as a controlling parameter, other than time, 
reaches its critical value. In our case, time is the controlling parameter and 
we are interested in the behavior of the solution as time reaches a critical 
value. 

The above arguments can be made more precise. The original system of 
equations for the modes F U G is given by 

*f- =R(t,«m 

where u = ({uk}), k G F U G is the vector of Fourier coefficients of u 
and R is the Fourier transform of the operator H. The system should be 
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supplemented with an initial condition u(0) = Uq. The vector of Fourier 
coefficients can be written as u = (u,u), where u are the resolved modes 
(those in F) and u the unresolved ones (those in G). Similarly, for the 
right hand sides (RHS) we have R(t,u) = (R(t,u), R(t,u)). Note that the 
RHS of the resolved modes involves both resolved and unresolved modes. In 
anticipation of the construction of a reduced model we can rewrite the RHS 
asR(t,u) = fl(°)(t,«) = (R^{t,u),R^(t,u)). For each mode u k , k G FUG, 
we can decompose R k (t,u) as 



Rf\t, U (t))=f2^ ) R ( a\t,u(t)). 



i=l 

Thus, the equation for the the mode u k , k € F U G is written as 

^ = R k (t,u) = R k °\t,u(t)) = f>f } Rf k \t, «(t)) (1) 

i=i 

Note that not all the coefficients af\ i = 1, . . . , m have to be nonzero. As is 
standard in renormalization theory [5] , one augments (with zero coefficients) 
the RHS of the equations in the full system by terms whose form is the 
same as the terms appearing in the RHS of the equations for the reduced 
model. Dimensional reduction transforms the vector = (a^f\ . . . ,am ) 

to a« = (4 1} ,... , o[j'). The reduced model for the mode ul, k € F is given 
by 

^ = Ri 1 \tMt)) = ±aVR$(t,m (2) 

i=i 

with initial condition u' k (0) = uq^- We emphasize that the functions , i = 
1, . . . , m, k G F, have the same form as the functions R^ , i = 1, . . . , m, k G 
F, but they depend only on the reduced set of modes F. This allows one to 
determine the relation of the full to the reduced system by focusing on the 
change of the vector to aS l \ Also, the vectors and do not have 
to be constant in time. This does not change the analysis that follows. 

Define m quantities Ei, i = l,...,m involving only modes in F. For 
example, these could be L p norms of the reduced set of modes. To proceed 
we require that these quantities' rates of change are the same when computed 
from (P| and P), i.e. 

dEi(u) dEi(u>) . 
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Note that similar conditions, albeit time-independent, lie at the heart of 
the renormalization group theory for equilibrium systems ([5] p. 154). In 
fact, it is these conditions that allow the definition and calculation of the 
(renormalization) matrix whose eigenvalues are used to calculate the critical 
exponents. In the current (time-dependent) setting, the renormalization 

matrix is defined by differentiating dEl ^ with respect to and using ([3]) 



to obtain 



d (dEi{u)\ 8 fdEi(u')\da { ^ 



daf\ dt J f^da^K dt J oaf)' 



i,j = 1, . . . ,m. (4) 



0a W 

We define the renormalization matrix M^j = — fey, k,j = l,...,m, as 

da) 

well as the matrices A^j = ~|oT ^ > k>J = 1j ■ ■ ■ ,m and Bkj = 

q %) d ^dt i " 1 ' = Equations (j3|) can be written in matrix 

form as 

A = MB (5) 

The entries of A describe the contributions of the different terms appearing 
on the RHS of the full system to the rate of change of E{ . The same for the 
entries of matrix B and the reduced model. 

The eigenvalues of the matrix M contain information about the behav- 
ior of the reduced system relative to the full system. In fact, they measure 
whether the full and reduced systems deviate or approach. In the renormal- 
ization theory of critical phenomena, the eigenvalues of M at the critical 
point are used to analyze the system properties close to criticality. The 
analysis is based on the assumption that the eigenvalues of M change slowly 
near the critical point so that even if one cannot compute exactly on the 
critical point, it is possible to get an accurate estimate of them by computa- 
tions near the critical point. Then, one performs a linear stability analysis 
near the fixed point and computes the system properties. The situation in 
the case of singularities of PDEs is different. In this case, the eigenvalues 
of M vary most rapidly near the singularity, due to the full system's rapid 
deterioration. Thus, we are not able to use linear stability analysis near the 
singularity. However, we are still able to extract the relevant blow-up rates 
(see Section [5|). 
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1.1 An instructive example 



We use the ID inviscid Burgers equation as an instructive example for the 
constructions presented in this section. The equation is given by 

ut + uu x = 0. (6) 

Equation ([6]) should be supplemented with an initial condition u(x, 0) = 
uq(x) and boundary conditions. We solve © in the interval [0, 2tt] with 
periodic boundary conditions. This allows us to expand the solution in 
Fourier series 

u M (x,t)= £ Mt)e ikx , 
keFuG 

where F U G = [— 4^ — 1]. We have written the set of Fourier modes 
as the union of two sets in anticipation of the construction of the reduced 
model comprising only of the modes in F = [-y,y — 1], where N < M. 
The equation of motion for the Fourier mode becomes 

du k ik m-^ 

-#=-2 L W ( 7 ) 

p+q=k 
p,q£FUG 



1.1.1 The t-model 

We need to choose a reduced model for the modes in F. We use a reduced 
model, known as the i-model, which follows correctly the behavior of the 
solution to the inviscid Burgers equation even after the formation of shocks 
[U [15] • The t-model was first derived in the context of statistical irre- 
versible mechanics |10| and was later analyzed in [U [15] . It is based on the 
assumption of the absence of time scale separation between the resolved and 
unresolved modes. We will use the same model for the case with nonzero 
viscosity and comment on its validity when appropriate. For a mode u' k in 
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F the model is given by 




2 u 'p u 'i 



p+q=k 
p£F, q&F 

ik 



p+q=k r+s=q 
p£F,q£G r<=F,seF 



- y Yl ~* 1" ^ U ' rU ' s U ' q ' ^ 



p+q=k L r+s=p 
p£G,q<=F reF,seF 



The first term on the RHS of ([8]) is of the same form as the first term in 
([7]), except that the term in (|8|) is defined only for the modes in F. The 
viscous term is the same. The third and fourth terms in ([8]) are not present 
in ([7j). They are cubic in the Fourier modes and they are effecting the drain 
of energy out of the modes in F. We should note here that the cubic terms 
in the t-model do not depend on the viscosity. To conform with the notation 
in Section Q] we rewrite as 




d 



ik 



p+q=k 
p£F,q£F 



a. 



2 



(1) 




ik 



p+q=k 
p£F, q£G 



r+s=q 
r£F,s£F 




Ik 



p+q=k 
p€G,q£F 



r+s=p 
r£F, seF 
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where = 1 and a)^ 1 = 1. We rewrite Equation ([7]) as 



,(!) 



UplLq 



p+q=k 
p,q£FuG 



+ 



2 ^ Up 

pG-FUG, gel 



2 ^ 



r+s=<j 
i-eFUG, sG-FUG 



p+q=k 
pGl,q&FUG 



2 ^ 



UfUg 



r+s=p 
rGFUG,sGFUG 



where of = 1 and = 0. The reader should note that we have introduced 
a new set of modes /. This is the set of unresolved modes for the full system. 
The reason for introducing the set / is that, as is the case in renormaliza- 
tion formulations, the terms appearing in the RHS of the equations at the 
different levels of resolution should be of the same functional form. The dif- 
ference between the different levels of resolution should be only in the range 
of modes used. Since the i-model involves a quadratic convolution sum with 
one index in the resolved range and the other in the unresolved range, we 
should use the same functional form when constructing the corresponding 
term for the full system. Thus, this term should involve a convolution sum 
with one index in the range F L) G and the other in /. 
Further, define 



ik 



p+q=k 
p,q£FUG 



and 



ik ^ — -v 
— > i 
2 ^ 

p+q=k 
pG-FUG, gGi" 



ik 
~2 



2 ^ 

Y. 

p+q=k 

pei,qeFuG 



r+s=q 
rGFUG, sGFUG 

ip 



-t- 



r+s=p 
rGFUG,sGFUG 



II,, 



9 



Also, define 
and 



ik 



p+q=k 
p,q€F 



ik \ - , 

y % 

p+q=k 
p£F,q£G 



J1 ' ' 

r+s=g 



p€G,q£F 



r+s=p 
r£F,s£F 



Thus, the equations of motion for the resolved modes in the full system and 
the reduced model can be written as 



and 



du k 
dt 



du 'k 
dt 



i=i 



i=i 



(9) 



(10) 



To proceed, we need to define the quantities Ei, i = 1, . . . , m. In our 
case, m = 2 and we need to define E\ and £2- The choice of the Ei is not 
unique. We chose for our experiments E± = ^ |itfc| 2 and £2 = Yl l^fcl 4 - 

The rates of change of the Ei are given for the full system by 



^ = YafhRe(R^(t,u(t))u* k ) + a^2Re(R^ (t,u(t))u* k ) 

k€F 



and 

dE 2 
dt 



Y^afhRe{2Rf k \t,u{t))\u k \^ 

k£F 



where u* k is the complex conjugate of u k - Similarly, for the reduced system 
we have 



^ = J2a^2Re(R^(t,u'(t))u' k *) + a^2Re(R^ (t,u'(t)K) 

k&F 
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and 



^ = ^aS 1) 2 J Re(2i?«(t,n'(t))|<| 2 n-) + a^2Re{2R^{t,u' (t))\u' k \^) 

The equations for the rates of change of the Ei can be used for the compu- 
tation of the 2x2 matrices A and B through the relations @ of Section 

m 



2 The mesh refinement algorithm 

We continue our presentation with the mesh refinement algorithm. The 
construction in the previous section requires the exact knowledge of an ac- 
curate reduced model. This means, the knowledge of both the functional 
form of the reduced model and the associated coefficient vector In fact, 
it is possible to relax this constraint by requiring the knowledge only of the 
functional form of the reduced model, i.e. knowledge of the vector but 
not of This can be considered as a time-dependent generalization of 
the Swendsen renormalization algorithm (e.g. see the nice presentation in 
Ch. 5 of [5]), even though here we do not have a statistical framework. The 
Swendsen algorithm is based on the observation that knowledge of only the 
functional form of the reduced model but not necessarily of the associated 
coefficient vector is enough for computing quantities of the reduced sys- 
tem. In particular, the matrix B can be calculated by using the resolved 
modes' values as computed from the full system. 

As we have mentioned before, the entries of B describe the contributions 
of the different terms appearing on the RHS of the reduced system to the rate 
of change of Ei (the same for the entries of matrix A and the full model) . The 
determinant of the matrix B measures whether there is need for the reduced 
system to transfer energy to smaller scales. The time instant when det B 
becomes nonzero, Tg, signals the onset of energy transfer from the modes in 
F to the modes in G. The determinant of the matrix A measures whether 
there is need for the full system to transfer energy to smaller scales. The time 
instant when det A becomes nonzero, Ta, signals the onset of underresolution 
of the full system. The time interval [Tb,Ta) is our window of opportunity to 
refine the mesh, without losing accuracy and without wasting computational 
resources. We will use the value of det B as a criterion to decide when it is 
time to refine the mesh. 

Note that if there exists a singularity, the interval AT = Ta — Tb will 
shrink to zero as we increase the resolution. The converse is not necessarily 
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true. If AT appears to converge to zero as we increase the resolution does 
not mean that there certainly exists a singularity. Since all the calculations 
are finite, there is only a maximum resolution that we can afford. It may 
well be that an even larger, and presently unattainable, resolution could 
reveal that there is no singularity. 

The mesh refinement algorithm is given by: 

Algorithm 

1. Choose a value for TOL. For this value of TOL run a mesh refine- 
ment calculation, starting, say, from N star t modes to Nfi na i modes. 
For example, let N s t ar t = 32 and double at each refinement until, say 
Nfinai = 256 modes. Record the values of the quantities E{, i = 
1, . . . , m when N = Nfi na i and \detB\ = TOL. Let's call this simula- 
tion 51. 

2. For the same value of TOL run a calculation with N s t ar t = Nfi na i 
modes (for the example N sta rt = Nfinai = 256). Record the values of 
the quantities Ei, i = 1, . . . , m when \detB\ = TOL. Let's call this 
simulation 52. 

3. Compare to within how many digits of accuracy the quantities Ei, i = 
1, . . . , m computed from 51 and 52 agree. If the agreement is to within 
a specified accuracy, say 5 digits, then choose this value of TOL. If 
the agreement is in fewer digits, then decrease TOL (more stringent 
criterion) and repeat until agreement is met. 

4. Use the above decided value of TOL to perform a mesh refinement 
calculation with a larger magnification ratio, i.e. a larger value for the 

ratio Nf ina i /N s t a rt ■ 

The agreement in digits of accuracy between 51 and 52 depends on the form 
of the terms chosen for the reduced model. Even though we do not know the 
coefficients of the reduced model, knowledge of the correct functional form of 
the terms can affect significantly the accuracy of the results. This situation 
is well known in the numerical study of critical exponents in equilibrium 
phase transitions (see e.g. Ch. 5 in [5]). 

2.1 How to compute the coefficients of the reduced model 

When we only know the functional form of the terms appearing in the re- 
duced model but not their coefficients it is not possible to evolve a reduced 
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system. We present a way of actually computing the coefficients of the re- 
duced model as needed. If the quantities Ei, i = 1, . . . . in are e.g. Lp norms 
of the Fourier modes, then we can multiply Equations ([2]) with appropriate 
quantities and combine with Equations (|3|) to get 

m 

1=1 

m 

i=i 

i=l 

where ujj , i,j = 1, . . . ,m are the new RHS functions that appear. Note 
that the RHS of the equations above does not involve primed quantities. 
The reason is that here the reduced quantities are computed by using the 
values of the resolved modes from the full system. The above system of 
equations is a linear system of equations for the vector of coefficients a™. 
In fact, the matrix of the system is the transpose B T of the matrix B. The 
linear system can be written as 

B T a^ = e (11) 

where e = {—^^--, ■ ■ ■ , ~~^~^)- This system of equations can provide us 
with the time evolution of the vector . 

The determination of coefficients for the reduced model through the 
system (jlip is a time-dependent version of the method of moments. We 
specify the coefficients of the reduced model so that the reduced model 
reproduces the rates of change of a finite number of moments of the solution. 
This construction can actually be used as an adaptive way of determining 
a reduced model if one has access to experimental values of the rates of 
change of a finite number of moments. Suppose that we are conducting a 
real world experiment where we can compute the values of a finite number 
of moments on a coarse grid only. Then we can use the system (|lip at 
predetermined instants to update a model defined on the coarse grid. Results 
of this construction will be presented elsewhere. 



dEi(u) 
dt 

dE 2 {u) 
dt 
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3 Numerical results for the inviscid Burgers equa- 
tion 



We present numerical results of the mesh refinement algorithm for the in- 
viscid Burgers equation with the initial condition u(x,0) = sin(x). This 
initial condition leads to a singularity forming at time T c = 1. Figure Q] con- 
tains results about the time spent between refinement steps and the time 
reached with the maximum allowed resolution. We start from a resolution 
Nstart = 32 and allow a maximum resolution of Nfi na i = 8192. We present 
results for two values of the tolerance TOL1 = 10 -16 and TOL2 = 10 -6 . 
When the tolerance criterion is less strict the algorithm can reach later times 
before running out of resolution. 

In Figure [2] we compare the velocity field produced by the algorithm 
with N start = 32, N fina i = 8192 and TOL1 = 1(T 16 with the velocity field 
produced by the algorithm with N star t = Nf ina i = 8192 and the same tol- 
erance. It is obvious that the results are in very good agreement. However, 
the mesh refinement calculation was about 240 times faster. The final time 
reached by the algorithm is T = 0.962. 

3.1 The direct approach to calculating the blow-up rate 

A mesh refinement algorithm can be used not only to approach a potential 
singularity but also estimate the rate at which the solution or some function 
of it blows-up. We restrict ourselves to the case of an algebraic (in time) 
singularity, meaning that some function of the solution diverges as ~ \T C — 
T|~ 7 , where 7 > 0. Let us assume for a moment that T c is known. One 
obvious way of estimating 7, is to run the mesh refinement algorithm and 
store the values of the blow-up quantity , say £ n , n = 1, . . . , N, and the 
instant T n at which each refinement took place. Then, one can plot (in log- 
log) the values of the blow-up quantity at the different refinement instants 
T n as a function of the distance from the singularity T c — T n and estimate 
the slope of the curve. That would provide us with the blow-up rate. Here 
we are interested in showing how the same estimate can be obtained using 
properties of a renormalization flow, i.e. a coarse-graining (dimensional 
reduction) procedure. Before we proceed, we have to address the issue of 
the value of T c which is, in general, unknown. Thus, the value of T c has to 
be calculated from the algorithm. It is simple to see that small errors in the 
estimation of T c can lead to huge errors in the estimation of the blow-up 
rate. One way of estimating T c is the following : for different choices of 
T c , plot, in log- log coordinates, the values of the blow-up quantity at the 
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Figure 1: (a) Time spent between refinement steps for different tolerance 
values, (b) Time reached with the maximum allowed resolution. 
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Figure 2: Comparison of the velocity field produced at the time of termina- 
tion of the mesh refinement algorithm for two different magnification ratios. 
The first simulation has N s t ar t = 32 and Nfi na i = 8192 while the second has 

N s tart = Nfinai = 8192. 

refinement instants T n as a function of the distance from the singularity 
T c — T n and pick the value of T c for which this plot is a straight line. This 
can be decided by monitoring the value of the correlation coefficient for a 
linear regression. 

We present results of the above construction for the inviscid Burgers 
equation with the initial condition u(x,0) = sin(x). This initial condition 
leads to a singularity forming at time T c = 1. Figure [3] shows the log-log plot 
of the maximum absolute value of the velocity gradient log (max jff | n ) and 
of the inverse distance from the singularity time (1 — T n ) _1 as recorded at the 
different refinement steps T n . The slope of the curve is 7 = 1 ± 1(T 8 . Note 
that the minute error in the estimate shows that the refinement algorithm 
keeps the calculation well-resolved even very close to the singularity. The 
calculations were performed using the mesh refinement algorithm of Section 
[2] with the refinement tolerance criterion TOL = det B set to 10 -10 . We 
should note that for this value of TOL, the value of det A for the full system is 
still much smaller than the double precision roundoff threshold of 10 -16 . For 
this calculation we set N s t ar t = 32 and Nfi na i = 131072 and the algorithm 
terminated at time T = 0.996. The mesh refinement is about 3000 times 
faster than a calculation with N s t ar t = Nfinai = 131072. 
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-log(l-Tn) 



Figure 3: Log-log plot of the maximum absolute value of the velocity gradi- 
ent max 1^:1^ and (1 T n ^ ^ for the different refinement steps (indexed by 
n). 

4 Numerical results for the supercritical focusing 
Schrodinger equation 

We continue with numerical results about the supercritical focusing Schrodinger 
equation. The focusing Schrodinger equation is given by 

i— + Au + \u\ 2a u = 0, where a > (12) 
ot 

The equation needs to be supplemented by an initial condition u(x, 0) = 
uo(x) and boundary conditions. It has been conjectured by Zakharov [22] 
that in d dimensions, when a > 2/d and for sufficiently large initial condi- 
tion, the solution of (|12p will blow-up at a finite time T and the behavior 
of the solution close to the blow-up time is given by 

u(x,t) = {(2k(T - i))-|^ +i ^Q((2K(T - t)r 1/2 |x|), 

where Q(£) is a complex- valued function with appropriate decay properties 
and k and uj are parameters to be determined. For the maximum of the 
solution we have 

max|n(x,t)| ~ (T - ty^ as t -> T. 

Although the mathematical theory is not yet complete, overwhelming evi- 
dence from numerical and formal analytical calculations suggests that the 
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conjecture is true. Here, we restrict attention to the ID case and to periodic 
boundary conditions in the domain [0, 2tt]. In the ID case, according to the 
conjecture the solution exhibits a algebraic finite time blow-up when a > 2. 
Here we present results for the case a = 3. In the numerical experiments we 
used the initial condition 

«o(x,0) = iAex.p(— (x — 7r) 2 ), 

for different values of A. For this initial condition we have max|ito(x)| = 
A at x = it. Figure H] shows the initial condition for A = 1.35 and the 
solution as computed by the mesh refinement algorithm with N s t ar t = 48 
and N f ina i = 10368. The tolerance criterion TOL = detB was set to 10~ 16 . 
The algorithm was implemented with the t-model for the reduced system as 
in the case of inviscid Burgers. 

Table Q] contains the estimated blow-up exponents for the maximum of 
the solution for different values of A. For A = 1.242 the mesh refinement 
algorithm does not run out of resolution which signals the absence of a 
singularity. For all the other cases and for N s t ar t = 48 and N^ na i = 10368, 
the mesh refinement algorithm was about 200 times faster than a calculation 
performed with N star t = ^ final = 10368. Unlike the case of the inviscid 
Burgers equation, here we cannot estimate beforehand the exact time T of 
the blow-up. We do that in the way proposed in the previous section. In 
particular, for different choices of T, we calculated the correlation coefficient 
of the linear fit (in log- log coordinates), of the values of the blow-up quantity 
as a function of the distance from the singularity T — t and picked the value 
of T for which the correlation coefficient is largest . For all the cases shown 
in Table Q] the correlation coefficient is about 0.999999999. The algorithm 
is able to approach the estimated singularity instant T to within 5 x 10~ 5 
units of time. The conjectured blow-up exponent for the maximum of the 
solution when a = 3 is l/2er = 1/6 ~ 0.1667. The relative deviation of the 
estimated values of the exponent relative to the conjectured value of the 
exponent is within 1% for all the values of A examined except for A = 1.8 
and A = 2. 

We would like to make a comment about the discrepancy for A = 1.8 
and A = 2.. It is to be expected that if one keeps the same maximum 
resolution while increasing the magnitude of the initial condition, i.e. the 
value of A, after some value of A the algorithm runs out of resolution before 
it can come close enough to the singularity for the asymptotic behavior to 
settle in. To elucidate this point we also ran the mesh refinement algorithm 
with N fi na i = 34992 for A = 1.8 and A = 2. The estimated values of the 
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Figure 4: Supercritical (a = 3) Schrodinger equation with max\uo(x)\ = 
1.35. 

blow-up exponent are included in Table Q] in parentheses. As we see, if one 
uses large enough resolution, the relative deviation of the estimated values 
of the exponent relative to the conjectured value of the exponent decreases 
again to within 1%. Note that for N star t = 48 and Nfi na i = 34992 the 
mesh refinement algorithm is about 400 times faster than a calculation with 
Nstart = N final = 34992. As expected, the acceleration factor increases when 
N final /N s t a rt increases. 

5 Calculation of the blow-up rate as a critical ex- 
ponent 

As we have said, we are also interested in showing how the blow-up rate 
estimate can be obtained using properties of a renormalization flow, i.e. a 
coarse-graining process. There are two ways to do that: i) Wilson's or "phase 
transition" approach, where one treats the singularity as a fixed point of a 
renormalization transformation and computes the blow-up rate by analysis 
in the vicinity of the fixed point, and ii) the Widom-Kadanoff or "scaling 
approach", where one assumes the existence of certain scaling laws in the 
vicinity of the singularity and then combines them to obtain the blow-up 
rate. 
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max un(x) 


Est. exp. a 


|Rel. dev.| (%) 


Sing, form 


1.242 








1.243 


0.1652 


0.90 


{T-t)- a 


1.250 


0.1648 


1.14 




1.255 


0.1654 


0.78 




1.260 


0.1649 


1.08 




1.300 


0.1655 


0.72 




1.350 


0.1662 


0.30 




1.500 


0.1678 


0.66 




1.600 


0.1691 


1.44 




1.800 


0.1727 (0.1684) 


3.60 (1.01) 




2.000 


0.1766 (0.1696) 


5.93 (1.74) 





Table 1: Estimated blow-up exponents for the supercritical (<r = 3) 
Schrodinger equation. The relative deviation is from the conjectured value 
of l/2fj = 1/6 ~ 0.1667. 

5.1 The "phase transition" approach 

The key idea is that a series of successive refinement steps (going to smaller 
and smaller scales) can be seen (approximately) as a coarse-graining process 
in reverse. Thus, one can run the mesh refinement algorithm, compute 
and store the coefficients of the reduced model at each refinement step and 
then use them to reconstruct the renormalization flow from smaller to larger 
scales. In this case, the smallest scale that the refinement algorithm reached 
is the starting scale of the renormalization flow. For the case of a time- 
dependent PDE the mesh refinement algorithm allows us to get closer and 
closer to the singularity instant T c . Thus, the renormalization procedure 
will take us further and further away from T c . 

There are two ways to show how the renormalization flow can be used 
to compute the blow-up rate. The first, the "phase transition" approach, 
assumes that the phase transition is a fixed point of the renormalization flow 
and proceeds with an analysis near the fixed point (e.g. pp. 124-127 in [5]). 
However, as we mentioned at the discussion after ([5]), we do not use a linear 
stability analysis because the eigenvalues of M vary most rapidly near the 
fixed point. Instead, we deal with the full (nonlinear) renormalization flow. 

The second way, the "scaling" approach, is just a manipulation of dif- 
ferent scaling laws assuming to hold asymptotically near the singularity. Of 
course, both lead to the same expression for the blow-up rate. We choose 
to present both since it elucidates further the connection between the tech- 
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-log(Ln) 

Figure 5: Log-log plot of the maximum absolute value of the velocity gradi- 
ent max |^r| n and. the inverse length scale of the reduced, system for the 
different renormalization steps (indexed by n). 

niques presented in this paper and those used in the theory of equilibrium 
phase transitions. 

We start our presentation of the blow-up rate calculation with the "phase 
transition" approach (see e.g. [5]). Let us suppose that near the singularity 
instant T c a quantity £ behaves as \T C — T| -7 . For the case of Burgers this 
would be the maximum of the velocity gradient, i.e. max|^^|. We want 
to find the value of 7. As we have said we assume that we have computed 
and stored a sequence of coefficients for the reduced model, the associated 
length scale, the value of the blow-up quantity and the time of occurrence 
of the refinement step. Then, by simply reversing the sequence indexing, we 
have the necessary quantities for the description of a renormalization flow 
which starts close to T c and moves further away with every coarse-graining 
step. Since every renormalization step brings us further away from the crit- 
ical point T c , the values of the blow-up quantity become smaller with every 
renormalization step. Thus, if we coarse-grain the length scale at which 
we probe the problem by a factor of b at each step (where b > 1), then 
£ n+ i = with /3 2 > 0. This implies £ n ~ Z^ 2 and thus /?2 can be com- 
puted from the refinement algorithm data collected. The coefficient of the 
reduced model which monitors the deviation of the full and reduced model 
will increase with each renormalization step, i.e. a n +i = otnb^ 1 , with f3\ > 0. 
This implies a n ~ In 1 and (5\ can also be computed from the collected data. 
Moreover, repeated application of the recursive relation for the coefficient 
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Figure 6: Log-log plot of the coefficient a 2 of the t-model and the length 
scale of the reduced system l n for the different renormalization steps (indexed 
by n). 




log(l-Tji) 



Figure 7: Log-log plot of the coefficient a 2 of the i-model and 1 — T n for 
the different renormalization steps (indexed by n). 
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a n gives a n = ao(6^ 1 ) n - This relation is the analog of the recursive relation 
derived in the theory of phase transitions by linearization of the renormal- 
ization flow around the critical (fixed) point. Here we did not resort to a 
linearization procedure. To proceed, we need to estimate the behavior of ao, 
the starting point of the renormalization flow. In the theory of phase tran- 
sitions, the behavior of the coefficient «o is assumed to be linear in \T C — T\. 
However, there is no a priori reason for such a behavior. We assume that 
otQ = C2\T C — T\ s , where 5 can also be computed from the collected data. 

Let us summarize what we have obtained so far. As we renormalize, the 
blow-up quantity decreases and the reduced model coefficient that monitors 
the deviation of the full and reduced model increases. Following the phase 
transition approach we thus assume that if we take enough renormalization 
steps then we have 



where u,v are quantities of the same order and C\,Ci are constants that 
depend on the initial conditions. We can eliminate n in the above two 
relations and get 



Thus, we have expressed the blow-up rate exponent 7 as a function of scaling 
exponents that are associated with properties of the renormalization flow. 

Before we conclude with this approach, we need to make one more com- 
ment. We have said before that the "phase transition" approach treats the 
singularity as a fixed point of the renormalization flow. To do that one has 
to construct a differential equation for the evolution of the coefficient a with 
respect to I. Note that by the way we have defined it, a is dimensionless. The 
equation for its evolution with changes in I is given by lOjr = (3(a) [5]. The 
RHS of the equation is called the beta function and its zeros determine the 
fixed points of the renormalization flow. Since a = for some constant 

C, we have l^f- = CftZ^ 1 . So, the beta function is 0(a) = Cp 1 ^ 1 = fta. 
So, the only fixed point of the beta function is a = 0. If ft > then a = 
corresponds to I = 0, i.e. the zero scale. But this is exactly the active scale 
reached at the instant that the singularity occurs. So, the singularity is in- 
deed a fixed point of the renormalization flow as long as ft > 0. Moreover, 
if ft > 0, this fixed point is unstable, so that if we start close to it, the 
renormalization flow will take us further away. This is indeed the case for 
the Burgers equation as we show numerically in the next section. 

This concludes the "phase transition" approach. 




u and C 2 \T C -T 



n 



= V 



£ ~ \T C — T\ 7 , with 7 = —5—. 
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5.2 The "scaling" approach 

We conclude with the "scaling" approach which is based on direct com- 
bination of the different scaling laws associated with the renormalization 
flow. Indeed, let £ ~ \T C — T| -7 , where 7' is the blow-up rate exponent 
to be estimated. If we assume that near T c we have £ ~ Z~^ 2 , a ~ l^ 1 
and a ~ \T C — T\ s , we can use the renormalization flow to estimate P\, P2 
and 5. Then a straightforward combination of the three scaling laws leads 
to 7' = j^-. So, 7' = 7 and as expected this approach leads to the same 
expression for the blow-up rate exponent as the "phase transition" approach. 

Figures [M3 show how one can use the above construction to estimate 
the blow-up rate 7 from renormalization flow quantities. Recall that the 
coefficient of the reduced model that monitors the deviation of the reduced 
and full systems is Also, that the index n appearing in the figures is 
used now to count the renormalization steps which are the opposite of the 
refinement steps. The length scale l n at which we probe the system for 
the different renormalization steps is the length scale of the reduced model. 
This means that if we have a full system calculation with N n modes, then 
l n = 2jt-, since the reduced model has half the resolution of the full system. 

From the data we estimate the exponents /?2 = 0.670±0.001, @\ = 0.739± 
0.007, and 6 = 1.1026 ± 10" 9 . From these estimates we get 7' = 1 ± 0.01. 
Thus, when we compute the blow-up rate using solely renormalization flow 
quantities, the estimation error is larger than when computing this rate 
directly. This is to be expected since we had to combine three empirically 
determined scaling laws, each one of which comes with its own error and 
also relies entirely on the adequacy of the reduced model. Nevertheless, the 
obtained accuracy is acceptable and moreover, it highlights the accuracy of 
the i-model for this equation. 

Finally, since (3\ = 0.739 > 0, we conclude that the singularity is an 
unstable fixed point of the renormalization flow (see discussion at the end 
of Section I5TT1) . 

6 Conclusions and future work 

We have presented a mesh refinement algorithm, inspired by renormalization 
constructions in critical phenomena, which allows the efficient location and 
approach of a possible singularity. The algorithm assumes knowledge of 
an accurate reduced model. In particular, it assumes knowledge of the 
functional form of the reduced model but not of the actual coefficients. We 
provide a way of computing the necessary coefficients on the fly as needed. 
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On a theoretical level, the algorithm can be used to study the behavior of 
(near-) singular solutions. On the practical side, it can be used as a mesh 
refinement tool. 

We have only examined the simple case of periodic boundary conditions 
and the mesh refinement performed was uniform. We plan to extend the con- 
structions presented here to a real space formulation which will allow the 
treatment of non-periodic boundary conditions and more complicated ge- 
ometries. In that case, one can divide the domain into sub-domains and ap- 
ply the mesh refinement algorithm individually in the different sub-domains. 
In addition, the algorithm can be modified to perform mesh-coarsening after 
the computationally intensive time interval of the simulation has passed. 

The original motivation behind the development of the algorithm was 
the open problem of the formation of singularities in finite time for the 
incompressible Euler and Navier-Stokes equations of fluid mechanics. In 
addition to helping with the issue of singularity formation, we hope that the 
algorithm can be of use in the simulation of real world flows by allowing a 
better assessment of the onset of underresolution. 
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